Analytic Delay Differential Equations

A delay differential equation is one whose derivatives are functions of the dependent variable at a different value of the independent variable. Typically they are integrated numerically by providing a history function of the same length as the delay. This presentation will consider equations that can be solved exactly without the need for any history function.

Consider first the equation

$f^{'} (t) = a f (t - T)$

where T is a constant displacement. The solution to the equation without displacement, $T = 0$ , is a simple exponential function. This function has the special feature of separability when its argument is the sum of two values. When solving the equation with displacement, the exponential in the independent variable will cancel from both sides leaving an exponential of the displacement on the right-hand side. The remaining transcendental equation is then soluble by the Lambert W function.

Assuming a solution of the form $e^{c t}$ , one has

$\begin{array}{l} c e^{c t} = a e^{c (t - T)} \\ c T e^{c T} = a T \end{array} \to c = \frac{W (a T)}{T}$

where the left-hand side has been arranged to match that of the defining equation of the Lambert function. The final solution to the differential equation is

$f (t) = exp [\frac{W (a T)}{T} t]$

which is easily verified. The branch of the Lambert function is left unspecified here, but if one wants to recover the expected solution as the delay goes to zero

$f (t) = e^{a t}, T \to 0$

then one can only use the principal branch, which is linear around the origin, since the other branches are all singular there. It should be noted than only $W_{0}$ and $W_{- 1}$ have portions of their domain where they remain real, so that in general this solution contains oscillatory behavior.

For $a = 1$ here is how the principal branch of the solution compares to an exponential:

This is just a simple scaling of the exponential because the combination $W (T) / T$ decreases monotonically from unity as T increases.

The solution generalizes easily to higher-order equations. Given the equation

$f^{(n)} (t) = a f (t - T)$

assume the same solution form and rearrange as before:

$\begin{array}{l} c^{n} e^{c t} = a e^{c (t - T)} \\ \frac{c T}{n} e^{c T / n} = \frac{T}{n} a^{1 / n} \end{array} \to c = \frac{n}{T} W (\frac{T}{n} a^{1 / n})$

The n independent solutions to this equation are given by introducing nth roots of unity multiplying the coefficient in the exponential,

$f_{k} (t) = exp [e^{2 π i k / n} \frac{n}{T} W (\frac{T}{n} a^{1 / n}) t], k = 1, 2 \dots, n$

which generally leads to oscillatory solutions for $k \neq n$ .

Now consider an equation with the constant displacement in the opposite direction:

$f^{'} (t) = a f (t + T)$

Following the procedure as before gives

$\begin{array}{l} c e^{c t} = a e^{c (t + T)} \\ - c T e^{- c T} = - a T \end{array} \to c = - \frac{W (- a T)}{T}$

with the solution

$f (t) = exp [- \frac{W (- a T)}{T} t]$

Since the principal branch of the Lambert function becomes complex when its function argument is less than $- e^{- 1}$ , this solution will become oscillatory for large enough delay. This is a marked difference from the first equation.

For $a = 1$ here is how the real (blue) and imaginary (purple) parts of the principal branch of the solution compare to an exponential (red):

The behavior changes significantly as already noted when $T > e^{- 1} \approx 0.368$ .

The solution again generalizes easily to higher-order equations. Given the equation

$f^{(n)} (t) = a f (t + T)$

again assume the same solution form and rearrange:

$\begin{array}{l} c^{n} e^{c t} = a e^{c (t + T)} \\ - \frac{c T}{n} e^{- c T / n} = - \frac{T}{n} a^{1 / n} \end{array} \to c = - \frac{n}{T} W (- \frac{T}{n} a^{1 / n})$

The n independent solutions to this equation are again given with nth roots of unity,

$f_{k} (t) = exp [- e^{2 π i k / n} \frac{n}{T} W (- \frac{T}{n} a^{1 / n}) t], k = 1, 2 \dots, n$

with the same general behavior.

Now consider an equation with equal delays in both directions:

$f^{'} (t) = a f (t - T) + b f (t + T)$

This equation again has an exponential solution, but the equation determining the coefficient in the exponential is no longer soluble with the Lambert function:

$c = a e^{- c T} + b e^{c T}$

One can however define a new double Lambert function as the solution to the nonlinear transcendental equation

$W (x, y) = x e^{- W (x, y)} + y e^{W (x, y)}$

The details of this function are considered in another presentation. The solution to the differential equation in terms of this double Lambert function is

$f (t) = exp [\frac{W (a T, b T)}{T} t]$

Here is how the real (blue) and imaginary (purple) parts of the principal branch of the solution compare to an exponential (red) for a variety of coefficient values:

This solution does not generalize simply to higher-order equations of the form

$f^{(n)} (t) = a f (t - T) + b f (t + T)$

because the equation determining the coefficient in the exponential is no longer soluble with the double Lambert function,

$c^{n} = a e^{- c T} + b e^{c T}$

since one cannot easily take an nth root of the linear combination of exponentials. Instead it appears one would have to define another new function that extends the double Lambert function.

The same appears to be true for an equation with unequal delays, such as

$f^{'} (t) = a f (t - T_{1}) + b f (t + T_{2})$

The equation determining the coefficient in the exponential is in this case

$c = a e^{- c T_{1}} + b e^{c T_{2}}$

which again is not soluble with the double Lambert function as defined. For any mixture of delay values and signs, one would have to adapt the numerical inversion used for the double Lambert function to each individual case.

To understand the limits of analytic solutions to delay equations, consider

$f^{'} (t) = a [f (t - T)]^{n}$

where n is an arbitrary value. The equation without delay is soluble by a binomial form:

$\begin{array}{l} f^{'} = a f^{n} \\ f^{'} f^{- n} = a \\ \frac{f^{1 - n}}{1 - n} = a t + C \end{array} \to f = [(1 - n) a t + f_{0}^{1 - n}]^{\frac{1}{1 - n}}$

It is easy to verify that this expression is a solution of the equation without delay. The expression represents a generalization of the exponential function to arbitrary powers of the right-hand side. From the limit definition of the exponential one has

$lim_{n \to 1} [(1 - n) a t + f_{0}^{1 - n}]^{\frac{1}{1 - n}} = lim_{n \to 1} f_{0} [1 + (1 - n) \frac{a t}{f_{0}^{1 - n}}]^{\frac{1}{1 - n}} = f_{0} e^{a t}$

for consistency. Unfortunately this expression does not have the special feature of separability of the exponential mentioned above, so a complete solution to the equation under consideration is not immediately available. Indeed, to see how complicated this equation is try a power series expansion to first order in T of its right-hand side:

$\begin{array}{c} f^{'} = a [f - T f^{'} + \dots]^{n} \approx a [f^{n} - n T f^{'} f^{n - 1}] \\ f^{'} f^{- n} + a n T \frac{f^{'}}{f} = a \\ \frac{f^{1 - n}}{1 - n} + a n T ln f = a t + C \\ \frac{f^{1 - n}}{a n T} + ln f^{1 - n} = \frac{1 - n}{n T} t + \frac{f_{0}^{1 - n}}{a n T} + ln f_{0}^{1 - n} \\ \frac{f^{1 - n}}{a n T} exp (\frac{f^{1 - n}}{a n T}) = \frac{f_{0}^{1 - n}}{a n T} exp (\frac{f_{0}^{1 - n}}{a n T}) exp (\frac{1 - n}{n T} t) \end{array}$

Now recognizing that the combination $z e^{z}$ is the inverse of the Lambert W function, one has

$f (t) \approx [a n T W [W^{- 1} (\frac{f_{0}^{1 - n}}{a n T}) exp (\frac{1 - n}{n T} t)]]^{\frac{1}{1 - n}}$

It should be emphasized that this is not an exact solution to the equation under consideration. The next approximation with a power series expansion to second order in T does not appear easily soluble. The binomial expression can be recovered from this one by setting $W (z) \approx ln z$ and $W^{- 1} (z) \approx e^{z}$ , although these approximations are not exact asymptotically.

Seeing this level of complexity from a first-order approximation gives one pause. Clearly some other analytic approach is needed.

One approach is Padé approximants, which are ratios of polynomials constructed to match the power series of a function to a given order. They are very often more accurate than the original power series, approximating the original function outside of the radius of convergence of the power series. This nature of stability makes them more useful for approximating unknown solutions to differential equations as well.

Begin with the simplest so-called diagonal approximant

$f (t) \approx \frac{f_{0} + a_{1} t}{1 + b_{1} t}$

For all equations in this presentation, the coefficient a represents a scaling of quantities with a dimension of time according to the order of differentiation in each equation. As such it can easily be restored if needed, so will be omitted in this section. Putting the approximant in the general delay differential equation under consideration, one has

$f^{'} (t) \approx \frac{a_{1} - b_{1} f_{0}}{(1 + b_{1} t)^{2}} \approx (\frac{f_{0} - a_{1} T + a_{1} t}{1 - b_{1} T + b_{1} t})^{n} \approx [f (t - T)]^{n}$

Since there are only two constants to determine, expand each side to first order in the temporal variable:

$(a_{1} - b_{1} f_{0}) (1 - 2 b_{1} t) \approx (\frac{f_{0} - a_{1} T}{1 - b_{1} T})^{n} (1 + \frac{n a_{1}}{f_{0} - a_{1} T} t - \frac{n b_{1}}{1 - b_{1} T} t)$

Equating coefficients gives the two equations

$a_{1} - b_{1} f_{0} = (\frac{f_{0} - a_{1} T}{1 - b_{1} T})^{n}$

$\frac{n b_{1}}{1 - b_{1} T} - 2 b_{1} = \frac{n a_{1}}{f_{0} - a_{1} T}$

While these equations are not particularly formidable, they cannot be solved exactly for an arbitrary power n, but one can naturally perform a numerical search for the root closest to the origin. A bit of numerical experimentation reveals that $b_{1}$ is always negative, so the search should start with that in mind. Taking $f_{0} = 1$ for convenience, the result is an approximation function that looks like this:

The function itself does not provide any sense of its accuracy as an approximation. A better indication is plotting the two sides of the differential equation above to see how well they match. With the left-hand derivative in red and the right-hand delayed function in purple, the result is

This first simple approximation is clearly of limited usefulness, but this is hardly unexpected.

The process can be extended with higher-order diagonal approximants for a more accurate and useful result. Consider the second-order diagonal approximant

$f (t) \approx \frac{f_{0} + a_{1} t + a_{2} t^{2}}{1 + b_{1} t + b_{2} t^{2}}$

Again set its derivative approximately equal to to the delayed power:

$\begin{array}{l} \frac{a_{1} - b_{1} f_{0} + 2 (a_{2} - b_{2} f_{0}) t + (a_{2} b_{1} - a_{1} b_{2}) t^{2}}{(1 + b_{1} t + b_{2} t^{2})^{2}} \\ \approx (\frac{f_{0} - a_{1} T + a_{2} T^{2} + (a_{1} - 2 a_{2} T) t + a_{2} t^{2}}{1 - b_{1} T + b_{2} T^{2} + (b_{1} - 2 b_{2} T) t + b_{2} t^{2}})^{n} \end{array}$

With the abbreviations

$\begin{array}{l} C_{0} = a_{1} - b_{1} f_{0} \\ C_{1} = \frac{2 (a_{2} - b_{2} f_{0})}{C_{0}} \\ C_{2} = \frac{a_{2} b_{1} - a_{1} b_{2}}{C_{0}} \end{array} \begin{array}{l} A_{0} = f_{0} - a_{1} T + a_{2} T^{2} \\ A_{1} = \frac{a_{1} - 2 a_{2} T}{A_{0}} \\ A_{2} = \frac{a_{2}}{A_{0}} \\ B_{0} = 1 - b_{1} T + b_{2} T^{2} \\ B_{1} = \frac{b_{1} - 2 b_{2} T}{B_{0}} \\ B_{2} = \frac{b_{2}}{B_{0}} \end{array}$

things are a bit more manageable:

$C_{0} \frac{1 + C_{1} t + C_{2} t^{2}}{(1 + b_{1} t + b_{2} t^{2})^{2}} \approx (\frac{A_{0}}{B_{0}})^{n} (\frac{1 + A_{1} t + A_{2} t^{2}}{1 + B_{1} t + B_{2} t^{2}})^{n}$

Setting $t = 0$ gives a nonlinear equation analogous to the first one in the previous approximation,

$C_{0} = (\frac{A_{0}}{B_{0}})^{n}$

which may then be canceled from the two sides of the equation. Since the Taylor series to the same order of two equal functions are necessarily identical, it is algebraically simpler at this point to rewrite the equation in the equivalent form

$(1 + C_{1} t + C_{2} t^{2}) (1 + B_{1} t + B_{2} t^{2})^{n} \approx (1 + b_{1} t + b_{2} t^{2})^{2} (1 + A_{1} t + A_{2} t^{2})^{n}$

In order to obtain three more equations to determine the four constants, expand each side to third order in the temporal variable. The two arbitrary powers have the same form under the rearrangement,

$\begin{array}{l} (1 + A_{1} t + A_{2} t^{2})^{n} \\ \approx 1 + n (A_{1} t + A_{2} t^{2}) + \frac{n (n - 1)}{2} (A_{1} t + A_{2} t^{2})^{2} \\ + \frac{n (n - 1) (n - 2)}{6} (A_{1} t + A_{2} t^{2})^{3} \\ \approx 1 + n A_{1} t + [n A_{2} + \frac{n (n - 1)}{2} A_{1}^{2}] t^{2} \\ + [n (n - 1) A_{1} A_{2} + \frac{n (n - 1) (n - 2)}{6} A_{1}^{3}] t^{3} \end{array}$

while the squared binomial to the same order is

$(1 + b_{1} t + b_{2} t^{2})^{2} \approx 1 + 2 b_{1} t + (b_{1}^{2} + 2 b_{2}) t^{2} + 2 b_{1} b_{2} t^{3}$

One merely need now equate coefficients of powers for the determining equations:

$n B_{1} + C_{1} = n A_{1} + 2 b_{1}$

$\begin{array}{l} [n B_{2} + \frac{n (n - 1)}{2} B_{1}^{2}] + n B_{1} C_{1} + C_{2} \\ = [n A_{2} + \frac{n (n - 1)}{2} A_{1}^{2}] + 2 n b_{1} A_{1} + b_{1}^{2} + 2 b_{2} \end{array}$

$\begin{array}{l} [n (n - 1) B_{1} B_{2} + \frac{n (n - 1) (n - 2)}{6} B_{1}^{3}] \\ + C_{1} [n B_{2} + \frac{n (n - 1)}{2} B_{1}^{2}] + n B_{1} C_{2} \\ = [n (n - 1) A_{1} A_{2} + \frac{n (n - 1) (n - 2)}{6} A_{1}^{3}] \\ + 2 b_{1} [n A_{2} + \frac{n (n - 1)}{2} A_{1}^{2}] + n A_{1} (b_{1}^{2} + 2 b_{2}) + 2 b_{1} b_{2} \end{array}$

Lovely! Clearly more complicated than the previous approximation, but still easily tractable numerically. The resulting function, again with $f_{0} = 1$ , looks similar to the previous case,

but again the function itself does not provide the information contained in the differential equation. With the left-hand derivative in red and the right-hand delayed function in purple, those two sides in comparison are

This is clearly a much better approximation than before, without an excess of added complexity. The approximation can be improved upon with a third- or fourth-order approximant, with a corresponding increase in algebraic complexity.

Padé approximants are really quite magical!

Uploaded 2021.04.10 — Updated 2022.12.17 analyticphysics.com